function [Y, periodStartConditions, spellDurations, spell2period] = generate_fullmodel_sequence_ctsTime2(indivConfigs, initialConditions, XbetaFixed, Xdynamic_staticInputs, paramsParts)
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	% This function generates a sequence of adoptions / exits according to a continuous-time model, and the timestamps of adoptions and
	% exits are stored.
	%
	%   - XbetaFixed is sums up the utility parts that do not vary at the (t-k) level (hence will not change
	%   as a function of laggedConsumers/ laggedProviders).
	%
	%  - Xdynamic_staticInputs contains the data parts at the (t-k) level (variables that do not depend on laggedConsumers and laggedProviders).
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	%%%%% Inputs:
	% indivConfigs:						object
	%	.consumerAdoptions:					object (as obtained by calling getConfig_consumerAdoptions)
	%	.providerAdoptions:					object (as obtained by calling getConfig_providerAdoptions)
	% initialConditions:				object:
	%	.laggedConsumersOrig:					K1 x 1
	%	.laggedConsumersDest:					K2 x 1
	%	.laggedProviders:						K x 1
	%	.popAtRisk:							object
	%		.providerAdoptions:					K x 1
	%		.consumerAdoptions:					K x 1
	%
	% XbetaFixed:						object:
	%	.consumerAdoptions:					K1 x K2 x T
	%	.providerAdoptions:					K x T
	%	.providerExits:						K x T
	%
	% Xdynamic_staticInputs:			object  (static inputs that, combined with laggedBases, serve to compute the Xparts
	%													that are impacted by lagged Y, --> these are the "static" things that do not change)
	%	.providerAdoptions:					object:
	%		.Xparts_static:						cell(NumParts, 1):
	%			{1}:								object:
	%				.X:									K x NumXstatic x T
	%				.X_FEs:								K x NumX_FEs x T
	%				.NumX_FE_vals:							1 x NumX_FEs
	%	.consumerAdoptions:					object:
	%		.Xparts_static:						cell(NumParts, 1):
	%			{2}:								object:
	%				.X:									K1 x NumXstatic x T
	%				.X_FEs:								K1 x NumX_FEs x T
	%				.NumX_FE_vals:							1 x NumX_FEs
	%			{3}:								object:
	%				.X:									K2 x NumXstatic x T
	%				.X_FEs:								K2 x NumX_FEs x T
	%				.NumX_FE_vals:							1 x NumX_FEs
	%
	% paramsParts:						object:
	%	.consumerAdoptions:					cell(NumParts,1)
	%	.providerAdoptions:					cell(NumParts,1)
	%	.providerExits:						cell(NumParts,1)
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	%%%%% Outputs:
	% Y:								object
	%	.consumerAdoptions:					sparse([NumSpells x (K1*K2)])
	%	.providerAdoptions:					NumSpells x K
	%	.providerExits:						NumSpells x K
	% periodStartConditions:			cell(NumPeriods,1):
	%	{tt}:								object:
	%		.laggedProviders:						K x 1
	%		.laggedConsumersOrig:					K x 1
	%		.laggedConsumersDest:					K x 1
	%		.popAtRisk:							obj:
	%			.providerAdoptions:					K x 1
	%			.consumerAdoptions:					K x 1
	%		.randomState:						obj (as obtained by calling rng)
	% spellDurations:                    NumSpells x 1
	% spell2period:                      NumSpells x 1  (gives value between 1 and NumPeriods)
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
		
	disp('Using the continuous-time model to generate a sequence with timestamps stored...');
	
	% Read dimensions
	[K,NumPeriods]       = size(XbetaFixed.providerAdoptions);
	[K1, K2, ~] = size(XbetaFixed.consumerAdoptions);
	
	%%%%% Initialize data structures to store results
	eventIdx = 0;
	eventTimeStamps = zeros(50000,1);
	eventTypes      = zeros(50000,1);
	eventDetails    = zeros(50000,3);
	
	%%%%% Initialize installed bases and populations at risk
	laggedConsumersOrig         = initialConditions.laggedConsumersOrig; % K1 x 1
	laggedConsumersDest         = initialConditions.laggedConsumersDest; % K2 x 1
	laggedProviders              = initialConditions.laggedProviders; % K x 1
	popAtRisk                 = initialConditions.popAtRisk; % object
	popAtRisk.providerExits      = laggedProviders; % K x 1
	
	%%%%% Loop over time periods
	for tt = 1:NumPeriods
		% Store conditions at start of period tt
		periodStartConditions{tt}.laggedProviders      = laggedProviders;
		periodStartConditions{tt}.laggedConsumersOrig = laggedConsumersOrig;
		periodStartConditions{tt}.laggedConsumersDest = laggedConsumersDest;
		periodStartConditions{tt}.popAtRisk         = popAtRisk;
		periodStartConditions{tt}.randomState       = rng;
		%%% Start day tt
		Tleft = 1;                           % time left until end of period tt
		
		while Tleft > 0
			% Get rates of arrival of events (updated every time a new event happens)
			installedBases.laggedConsumersOrig = laggedConsumersOrig;
			installedBases.laggedConsumersDest = laggedConsumersDest;
			installedBases.laggedProviders      = laggedProviders;
			arrivalRates_tt = getArrivalRates(indivConfigs, popAtRisk, paramsParts, XbetaFixed, Xdynamic_staticInputs, installedBases, tt);
			
			% Get sum of arrival rates (across all types of events)
			lambda_consumerAdoptions = arrivalRates_tt.consumerAdoptions;     % K1 x K2
			lambda_providerExits      = arrivalRates_tt.providerExits;           % K  x 1
			lambda_providerAdoptions  = arrivalRates_tt.providerAdoptions;       % K  x 1
			
			lambda_ttl_consumerAdoptions = sum(lambda_consumerAdoptions(:)); % 1 x 1
			lambda_ttl_providerExits      = sum(lambda_providerExits);      % 1 x 1
			lambda_ttl_providerAdoptions  = sum(lambda_providerAdoptions);  % 1 x 1
			
			lambda_total = lambda_ttl_consumerAdoptions + lambda_ttl_providerAdoptions + lambda_ttl_providerExits; % 1 x 1
			
			% Simulate time to next event
			A = draw_Exponential(lambda_total, 1); % 1 x 1
			
			% If event happens before end of the day: record it, update everything, reduce time left before end of day and iterate
			if A < Tleft
				evt_date = tt - Tleft + A;
			
				% Figure out what type of event happened: provider adoption, provider exit or consumer adoption (multinomial probs)
				cond_probs = [lambda_ttl_consumerAdoptions; lambda_ttl_providerAdoptions; lambda_ttl_providerExits]./lambda_total; % 3 x 1
				evt_type = find(mnrnd(1,cond_probs)==1); % integer between 1 and 3
				clear cond_probs;
				
				%%% Case of consumerAdoption
				if evt_type == 1
					% Draw origin of event
					cond_probs = sum(lambda_consumerAdoptions,2)./lambda_ttl_consumerAdoptions; % K1 x 1
					k1 = find(mnrnd(1,cond_probs)==1); % index between 1 and K1
					clear cond_probs;
					% Draw destination of event
					cond_probs = lambda_consumerAdoptions(k1,:)'; cond_probs = cond_probs./sum(cond_probs); % K2 x 1
					k2 = find(mnrnd(1,cond_probs)==1); % index between 1 and K2
					clear cond_probs;
					
					% Update installed base and popAtRisk
					laggedConsumersOrig(k1)             = laggedConsumersOrig(k1) + 1;
					laggedConsumersDest(k2)             = laggedConsumersDest(k2) + 1;
					popAtRisk.consumerAdoptions(k1) = popAtRisk.consumerAdoptions(k1) - 1;
				end
				
				%%% Case of providerAdoption
				if evt_type == 2
					% Draw location of event
					cond_probs = lambda_providerAdoptions./lambda_ttl_providerAdoptions; % K x 1
					kk = find(mnrnd(1,cond_probs)==1); % index between 1 and K
					clear cond_probs;
					
					% Update installed base and popAtRisk
					laggedProviders(kk)             = laggedProviders(kk) + 1;
					popAtRisk.providerAdoptions(kk) = popAtRisk.providerAdoptions(kk) - 1;
					popAtRisk.providerExits(kk)     = popAtRisk.providerExits(kk) + 1;
				end
				
				%%% Case of providerExit
				if evt_type == 3
					% Draw location of event
					cond_probs = lambda_providerExits./lambda_ttl_providerExits; % K x 1
					kk = find(mnrnd(1,cond_probs)==1); % index between 1 and K
					clear cond_probs;
					
					% Update installed base and popAtRisk
					laggedProviders(kk)             = laggedProviders(kk) - 1;
					popAtRisk.providerAdoptions(kk) = popAtRisk.providerAdoptions(kk) + 1;
					popAtRisk.providerExits(kk)     = popAtRisk.providerExits(kk) - 1;
				end

				% Reduce time left until end of day
				Tleft = Tleft - A;
			
			% If event happends after end of day: discard it, declare day over
			else
				evt_date = tt;
				evt_type = 0; % Day change
				
				% Reduce time left until end of day
				Tleft = 0;
			end
			
			% Store timestamp, type and details of event
			eventIdx = eventIdx + 1;
			if eventIdx > size(eventTimeStamps,1)
				eventTimeStamps = [eventTimeStamps; zeros(50000,1)];
				eventTypes      = [eventTypes;      zeros(50000,1)];
				eventDetails    = [eventDetails;    zeros(50000,3)];
			end
			eventTimeStamps(eventIdx) = evt_date;
			eventTypes(eventIdx)      = evt_type;
			if evt_type == 1   % Consumer adoption
				eventDetails(eventIdx,2) = k1;
				eventDetails(eventIdx,3) = k2;
			else if evt_type == 2 || evt_type == 3 % Provider adoption or provider exit
				eventDetails(eventIdx,1) = kk;
			end; end;
					
			% Clean up
			clear evt_date evt_type k1 k2 kk;
		end
		displayRemainingTime(tt, NumPeriods, 20);
	end
	
	% Memory pre-allocated: only keep those events that occurred (remove extra rows allocated)
	eventTimeStamps = eventTimeStamps(1:eventIdx);
	
	% Check that events appear in increasing order
	assert(all(eventTimeStamps(2:end) > eventTimeStamps(1:end-1)));
	
	% Compute spellDurations and spell2period
	spellStarts = [0; eventTimeStamps(1:end-1)];
	spellEnds = eventTimeStamps;	
	spellDurations = eventTimeStamps - spellStarts;        % NumSpells x 1
	spell2period = make_integers(floor(spellStarts) +  1); % NumSpells x 1 (gives values between 1 and NumPeriods)
		
	% Put together Y.consumerAdoptions, Y.providerAdoptions and Y.providerExits
	NumSpells = length(eventTimeStamps);
	Y.consumerAdoptions = sparse(NumSpells, K1*K2);
	Y.providerAdoptions  = sparse(NumSpells, K);
	Y.providerExits      = sparse(NumSpells, K);
	for uu = 1:NumSpells
		if eventTypes(uu) == 1 % consumer adoption
			k1 = eventDetails(uu,2);
			k2  = eventDetails(uu,3);
			myidx = sub2ind([K,K],k1,k2);
			Y.consumerAdoptions(uu, myidx) = 1;
		end
		if eventTypes(uu) == 2 % provider creation
			kk = eventDetails(uu,1);
			Y.providerAdoptions(uu,kk) = 1;
		end
		if eventTypes(uu) == 3 % provider exit
			kk = eventDetails(uu,1);
			Y.providerExits(uu,kk) = 1;
		end
		clear k1 k2 kk myidx;
	end
end
